Space-charge mechanism of aging in ferroelectrics: an exactly solvable 

two-dimensional model 
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A mechanism of point defect migration triggered by local depolarization fields is shown to explain 
some still inexplicable features of aging in acceptor doped ferroelectrics. A drift-diffusion model 
of the coupled charged defect transport and electrostatic field relaxation within a two-dimensional 
domain configuration is treated numerically and analytically. Numerical results are given for the 
emerging internal bias field of about 1 kV/mm which levels off at dopant concentrations well be- 
low 1 mol%; the fact, long ago known experimentally but still not explained. For higher defect 
concentrations a closed solution of the model equations in the drift approximation as well as an 
explicit formula for the internal bias field is derived revealing the plausible time, temperature and 
concentration dependencies of aging. The results are compared to those due to the mechanism of 
orientational reordering of defect dipoles. 

PACS numbers: 77.80.Dj,77.80.Fm,77.84.Dy,61.72.jd 
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I. INTRODUCTION 

The phenomenon of gradual change of physical prop- 
erties with time, called aging, is for a long time known 
feature of ferroelectrics, especiall y w hen acceptor doped 
EE 0, 0, S !, IS 0, 0, S El HO, Hi 0- Aging reveals 
itself in quasi-logarithmic decrease of the dielectric con- 
stant with time [l], 0] , reduction of domain wall mobility 
leading to stabilization in the aged domain structure [2J , 
altered shape of polarization loops in bothpoled and un- 
poled aged samples 0, 0, 0, US HH, HI, and related 
indications. A characteristic aging time, r, a clamping 
pressure on the domain walls, P c i, and an internal bias 
field, i?if,, were introduced as parameters quantifying ag- 
ing&Q. 

In past three decades several concepts were developed 
[3 fig to explain the aging phenomena in terms of do- 
main splitting [2], space charge formation [3, 14, elec- 
tronic charge tra ppin g at domain boundaries [9lll0|. ionic 
la. \YA [l8j or reorientation of defect dipoles 
12j. The latter concept based on the widely 



drift 
0-0. 

recognized mechanism of gradual orientation of defect 
dipoles formed by the charged acceptor defects and oxy- 
gen vacancies has suggested probably most successful 
quantitative explanation of many features relevant to ag- 
ing and fatigue in ferroelectrics, particularly, plausible 
time and temperature dependencies of En,. Neverthe- 
less, some long standing questions remain open, most 
pronounced among them the dependence of r and En, on 
the defect concentration d 0]. Resulting from the in- 
dividual cage motion of an oxygen vacancy around an 
acceptor defect insensitive to presence of other defect 
dipoles in the dipole reorientation model @, 0] r is ex- 
pected to be independent on the acceptor concentration. 
Similarly, independent contributions of different defect 



dipoles to the clamping pressure in this model have to 
result in En, directly proportional to the concentration. 
However, the experimentally observed saturation of En, 
at medium concentrations < 1 mol% as well as distinct 
concentration dependence of r [1 0, 0] provides indica- 
tions of some collective mechanism of aging. 

In this paper, therefore, we prove an alternative, space 
charge mechanism to explain the above mentioned fea- 
tures of aging. It may also be related to the self- 
polarization phenomenon and internal field-induced, mi- 
gratory polarization observed in thin ferroelectric films 
[lOL . Recently, a model quantifying the space charge 
mechanism was advanced (20 . l2ll . [22j which shows that 
the clamping pressure P c i ~ 1 MPa and the field En, ~ 
1 kV/mm comparable with experiments can result from 
the formation of space charge zones near charged do- 
main boundaries assuming small but finite mobility of 
charged defects. Aging was studied first for low defect 



, [21| , and the effect 
Here we apply the 



concentrations about 0.01 mol% 
of anisotropy was considered [22 1. 
isotropic, two-dimensional version of this model [2l| to 
study aging of unpoled ferroelectrics for a wide range of 
acceptor concentrations and present numerical and an- 
alytical results on temperature, concentration and time 
dependencies of En,. Gradual change of material proper- 
ties under the effect of the external dc or ac electric field, 
e.g. fatigue phenomenon, is not considered at this stage. 



II. 



MODEL OF A FERROELECTRIC GRAIN 



'Electronic address: yugenen@tgm.tu-darmstadt.de 



Main assumptions of the model [2l| used here are: (a) 
availability of mobile charged defects in amount sufficient 
to substantially compensate the bound charge at the do- 
main boundaries and (b) presence of strong local depo- 
larization fields in the unpoled ferroelectric material. 

Let us consider the assumption (a). The conductivity 
of perovskites has been extensivel y s tudied during the 
past two decades [H, Hi 111 US ESIH IH ■ It was estab- 
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lished that, depending on temperature and partial pres- 
sure of oxygen, perovskites may exhibit ionic or electronic 
conductivity, so that the n-type conductivity prevails un- 
der reducing conditions and the p-type in oxidizing atmo- 
spheres [23], [2(| [13] ■ The ionic conductivity prevails be- 
tween the above two areas of electronic conductivity and 
is clearly dominated by oxygen vacancies. It was shown 
also that the bulk mobilities of electrons and holes are not 
activated in BaTi03 and exceed the activated mobility 
of vacancies by many orders of the magnitude. Neverthe- 
less, at atmospheric partial oxygen pressure and not very 
high temperatures, the concentration of electronic carri- 
ers remains less than the concentration of vacancies by 
many orders of the magnitude and is by far not sufficient 
to compensate the polarization bound charge. Moreover, 
the bulk concentration of electronic carriers is so small 
that the corresponding Debye screening length strongly 
exceeds the typical grain size. This allows one to entirely 
neglect electronic screening of the local bound charges. 
The same applies also to ambient electronic carriers in 
the intergranular space though their concentration may 
exceed the bulk one by few orders of the magnitude [28| • 
The balance between the electronic and ionic species may 
be substantially changed locally right near the charged 
domain boundaries because of the electronic band bend- 
ing by very strong local depolarization fields; however, 
this possibility depends on the reduction potentials of 
certain metal dopants [301 ] . In our study we suppose for 
simplicity that concentration of electronic carriers is neg- 
ligible in and outside the ferroelectric grains. The oxygen 
vacancy concentration is usually fixed by acceptor defects 
even in the nominally undoped materials, since these are 
as a rule unintentionally acceptor doped [26|, [13, HI] . For 
these reasons, we will assume in the following that only 
oxygen vacancies participate in charge migration. 

Consider now the assumption (b) . In the perfectly or- 
dered domain system the bound charges at the domain 
boundaries can be fully compensated resulting in com- 
plete suppression of the depolarization fields inside the 
bulk ferroelectric. Experimental studies show that it is 
not the case in real, disordered systems. Considering the 
high resolution images of domain patterns in BaTi03 and 
lead zirconate titanate (PZT) samples by scanning trans- 
mission electron microsco py 13 II) , by bright field transmis- 
sion electron microscopy [321 ] and by secondary electron 
spectroscopy |33j one can observe numerous places where 
the uncompensated bound charges should emerge and, 
consequently, local depolarization fields can be present. 
Three relevant circumstances can be identified, namely, 
when a domain array a) ends up in an unpolarized area 
inside the grain, b) meets a domain wall of another do- 
main array instead of charged domain faces, c) ends up at 
the grain boundary contacting the wide enough unpolar- 
ized integranular area. To investigate quantitatively the 
consequencies of the presence of uncompensated bound 
charges in the material the third of the above mentioned 
cases, (c), will be representatively studied in the follow- 
ing. Since the depolarization field of a domain array de- 



cays exponentially at the distance equal to the domain 
width a [13, HU , it suffices to assume the intergrain spac- 
ing larger than a to get the neighbour grains electrically 
independent. That is why, for studying the effects of lo- 
cal depolarization fields, a model of a single ferroelectric 
grain surrounded by unpolarized medium may be chosen. 

Let us start with a quadratic ferroelectric grain of zero 
total polarization surrounded by an infinite dielectric 
medium. We imagine a two-dimensional periodic array of 
domains in the grain cut by the interfaces with dielectric, 
z = and z = L, perpendicular to the direction of spon- 
taneous polarization which is along the z-direction of a 
Cartesian coordinate system x,y,z (Fig. [TJ) . The system 
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FIG. 1: Scheme of a 2D-array of 180°-domain walls in 
a quadratic ferroelectric grain occupying the space | as | < 
L/2, < z < L. Straight arrows show the direction of the 
polarization and curved arrows the approximate direction of 
the local electric field. 

is supposed to be uniform in the y-direction so that no 
variable is y-dependent. If the length of the domains L 
along the z-axis is much larger than their width a along 
the x-axis, which is typically the case in experiments, 
electric field lines are effectively closed at the same side 
of the grain. Finite-element simulations of the electric 
field in the finite domain array shows a virtually periodic 
field pattern along the x axis everywhere but the very 
ends of the array [2(| HH . That is why we can consider a 
periodic domain array infinite along the x axis as a rep- 
resentative model for a finite multi-domain grain. This 
model configuration is well-known in the physics of po- 
larized media and was used for the study of equilibrium 
and dynamic properties of ferromagnetic [34 , |35| and fer- 
roelectric [36| materials. 

Furthermore, since both components of the depolar- 
ization field exponentially decrease towards the interior 
of the grain along the z-axis [13, HH, transport of the 
charged defects driven by the field is expected to occur 
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in the vicinity of the grain boundaries z = and z = L. 
Considering charge migration near the interface z — 
we can therefore assume the domains to be infinite along 
the z-axis without introducing a substantial error. Ap- 
parently, the same process of charge separation occurs at 
the other grain boundary, z — L, too. Calculating the 
forces exerted upon the domain walls, both ends of the 
domains must be taken into account. 

Hence, to study the depolarization field induced charge 
migration it is sufficient to consider the interface z = 
between the domain array occupying the ferroelectric 
half space z > and the dielectric medium occupying 
the half space z < 0. For simplicity, both media are 
assumed to be isotropic and characterized by the relative 
dielectric permittivities £/ and Ed, respectively. Due to 
polarization, the domain faces at z = are alternatively 
charged with the bound surface charge density [3|| 



p b (x,z) = a p 5(z) ]T (-1)" 

Tl— — OO 

(| + an-x^j , 



■ (-> ( — — an + x 



(1) 



where a p — 



P s | is the spontaneous polarization value, 
5(z) and 9{x) are the Dirac 8- function and the Heaviside 
unit step function, respectively. The choice of the origin 
in the middle of the positively charged domain face makes 
the problem also bilaterally symmetrical. 

Migration of charge carriers is governed by the drift- 
diffusion equation which is often used for description of 
electronic transport in semiconductors (37j : 



d t c= -V(/icE) + DAc, 



(2) 



where E(x, z, t) is the local electric field at the time t, 
c(x,z,t) is the concentration of positively charged, mo- 
bile species, /i and D are their mobility and diffusivity, 
respectively. We assume, for simplicity, that the latter 
two quantities are isotropic and connected by the Ein- 
stein relation, D = fjtkT '/ 'qf with k the Boltzmann con- 
stant, T the absolute temperature and qf the charge of 
the carriers. 

The electric field E(x,z,t) is determined self- 
consistently by the charged faces of the domains and 
the imbalanced charge density in the bulk p(x, z,t) = 
q/ [c(x, z, t) — cq] through Gauss' law 
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(3) 



where Co is the background concentration of the immobile 
acceptor defects warranting total electroneutrality, and 
£o is the permittivity of vacuum. 

For the boundary conditions to the system of equations 
([2]) and ([3]) we assume usual boundary conditions for the 
electric field at the interface between the two media [35[ 
and chemical isolation of the grain determined by the 
requirement of vanishing particle current at z = 0: 



Note that the latter condition is not crucial for the 
present model. Transport of charged species in the di- 
electric medium (intergranular space) and also through 
the grain boundary may be additionally included. Here, 
for simplicity, we consider only redistribution of mobile 
defects inside the grain. 

In the initial state, the system is locally neutral as- 
suming c(x, z, 0) = Co, while the electric field E(x, z, 0) = 
E°(cr p |x,z) is determined at the beginning solely by the 
surface charge density (fT|) and reads, as was calculated 
in Ref. (H, 



ne Q (e f + e d ) 
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cosh(7rz/a) + sin(7ra/a) 
cosh(-7rz/a) — sin(7ra;/a) 
I f cos(7rx/a) 



7T£ (£/ + Ed) 



■ arctan 



sinh(7rz/a) 



(5) 



fxcE z - Dd z c = . 



(4) 



The latter equations together with boundary conditions 
and Eqs. and (O complete the formulation of the 
problem which can now be treated numerically 



III. NUMERICAL EVALUATION OF THE 
INTERNAL BIAS FIELD 

Numerical solution of the problem is performed using a 
simple direct Euler scheme described in Ref. [2l| . Due to 
the periodicity and the bilateral symmetry of the task, it 
is sufficient to consider the charge redistribution within 
the area < x < a contaning a single domain wall at 
the plane x = a/2. Space and time are discretized. At 
every time step, the change in the carrier concentration is 
computed from the previous values of the concentration 
and the electric field using Eq. ((2]). Updated values of the 
field are then calculated using Eq. (J3J). The evaluation is 
repeated until convergence is achieved. 

As an example, we consider now the aging process in 
unpoled BaTi03 doped with a bivalent acceptor, e.g. Ni, 
Ca or Mg. For numerical simulations, the material pa- 
rameters at temperature 45° C are taken from Refs. [7l.[38j 
with the aim to compare results directly with those of 
the theory of dipole reorientation. Namely, it is assumed 
that P s = 2.7 • KT 1 C/m 2 , e f = 200, a = 0.5 ^m. Pos- 
itively charged oxygen vacancies are assumed as mobile 
species with qf twice the elementary charge and the mo- 
bility fi = 8.4 • 10 -22 m 2 / Vs implying activation energy 
of E a = I. leV @,m,[22l,[2p. 

Note that, for the dielectric permittivity the intrinsic, 
lattice value is taken [38| since the charge migration oc- 
curs at the mesoscopic scale of the domain width. The di- 
rectly measurable macroscopic permittivity may achieve 
much higher values due to the large contribution of the 
domain walls [39j which does not apply in the considered 
problem. Account of anisotropy of crystalline BaTi03 
would result in enhancement of the field magnitude by 
the factor ^JeajEc and in reduction of the field penetra- 
tion depth by the same factor [2l|, |22j, where e a — 2180 
and e c = 56 are the principal values of the permittivity 
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tensor [4(|. Since the electrostatic energy is quadratic in 
the field this would entail the increase in the values of the 
cl ampin g pressure and the internal bias field by the factor 
\/e a /e c ~ 6. However, for simplicity and comparability 
with the results of the dipole reorientation model @, [1] , 
we assume here the above introduced isotropic permittiv- 
ity e f . For the dielectric medium outside the ferroelectric 
grain we take the same but non-polarized material with 
£<* = £/. 

The system reveals two typical time scales: the drift 
time t m = a/fxE® ~ 7.8 • 10 6 s and the diffusion time tjj = 
a 2 /D ~ 2.2-10 10 s so that the ratio = t^/td ~ 3.6-1CT 4 
characterizes the contribution of diffusion to Eq. ([2]) [2l[ . 
Though very small, this contribution cannot be neglected 
because it influences the structure of space charge zones 
and provides compatibility with the boundary condition 
((4]). However, because of very small /3, large gradients 
of concentration arise near the negatively charged face of 
the domain at z = 0, a/2 < x < a which make the numer- 
ical procedure unstable. Since the parameter (3 affects 
nothing but the thickness of the positively charged layer 
piled up near the negatively charged domain face [2lj we 
assume in computations /3 = 5 ■ 10~ 2 keeping in mind 
that this may lead to mistakes if the thickness of the 
space charge zone in front of the positively charged do- 
main face becomes comparable with (3a. Details on the 
field and concentration profiles and their evolution with 
time are exemplary presented in Ref. [2l| for low dopant 
concentrations about cq — 0.01 mol%. Here computa- 
tions are extended over the region from cq — 0.01 mol% 
to 1 mol%. 

Having the charge density and the electric field calcu- 
lated, the time dependent forces exerted upon domain 
walls can be evaluated. The loss of domain wall mo- 
bility, characteristic of aging, results from relaxation of 
the energy of the electrostatic depolarization field due to 
piling up of the charged defects at the charged domain 
faces. Distribution of the energy of this field and of the 
consequent clamping pressure along the domain wall is 
very nonuniform, peaking at the domain boundaries [2lj . 
The average clamping pressure preventing the displace- 
ment of the domain wall from the energy minimum and 
the corresponding internal bias field may be estimated as 
follows. 

The thermodynamic force exerted upon the domain 
wall can be defined as the derivative of the energy of the 
electrostatic field on the domain wall diplacement. This 
derivative may be roughly estimated from the difference 
between the energy of the initial state of the system dis- 
played in Fig. [1] and the fully polarized state achieved 
by the virtual displacement of domain walls over the dis- 
tance of a/2. 

The energy of the electrostatic field per one half of the 
domain length reads 



1 r L / 2 r a 
W[E] = -e e f / dz dxE 

2 J-oc Jo 



(6) 

-oo Jo 

while the other half of the domain at z > L/2 contributes 



the same amount of the energy for symmetry reasons. In 
the virgin state of the grain, the depolarization field is 
represented by Eqs. j5|), and the energy of this field per 
one period of the structure is given by the well known 
formula Hi, 11 



W[E°(a p )} =0.85 



2 2 
47T£ £ / ' 



(7) 



In the course of aging, the electric field transforms to 
E = E°(er p ) + AE, where AE is the contribution to the 
field due to the charge carrier migration. The aged state 
of the structure in Fig. [T] serves as the initial state with 
the energy W[E] in the clamping force calculation. 

Consider now the virtual pairwise displacement of the 
domain walls to each other until they meet which leads 
to the full polarization of the system in the positive z 
direction. The resulting uniform bound charge at the 
grain boundaries at z = and z = L generates then the 
uniform depolarization field E^ = (0, 0, — P s /eo£f), < 
2 < L, so that the total electric field becomes equal to 
Ey = Erf + AE. The energy of the electrostatic field 
in this state amounts to W[E/]. The clamping pres- 
sure on the domain walls related to aging is provided by 
the time-dependent part of the energy difference VF[E/]- 
W[E], namely 



AW C , 



s s f 



L/2 



dz 



dxfEj-E ) AE 



(8) 



Calculating here the energy gain near the domain bound- 
ary at z = the integration limit L/2 can be extended 
to infinity because of the exponentially fast descrease of 
the fields E° and AE in both directions of z axis. 

Finally, the internal bias field in the z direction can be 
evaluated comparing the force 2P S £L [4l[ exerted upon 
the domain wall by an external field £ with the clamp- 
ing force 2AW c i/(a/2) accounting now for both domain 
boundaries. This results in the estimation 



E,< 



aLP, 



(9) 



Evaluation of the time-dependent field En, assuming 
the typical length of the domain wall L = 20a is shown 
in Fig. [2] for different dopant concentrations. All the 
curves demonstrate saturation of En, after a some charac- 
teristic (aging) time r which decreases with the growing 
concentration. The final, equilibrium value of En, first in- 
creases with the concentration but then levels off well be- 
low Co — 1.0mol% as is seen in Fig.[3J The inversed aging 
time t -1 rises almost linearly with the increasing doping 
in the whole range of concentrations studied. Saturating 
dependence of En, on Co is clearly seen in experiments on 
the acceptor doped PZT 0,0]; indications of both above 
described En, and r dependencies are observable on the 
acceptor doped barium titanate too 

Saturation of E^ as well as decreasing of aging time 
with growing concentration cq has a simple physical rea- 
son. The depth A of the space charge area emerging 
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FIG. 2: Eib as a function of time for variable acceptor con- 
centration Co which increases monotonically from the lower 
to the upper curve (solid lines). The dashed line presents the 
theoretical limit of high cq (see text). 
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FIG. 3: Saturated value of En and inverse aging time, r _1 , 
determined from the inflection point of plots in Fig. [2] are 
presented as functions of acceptor concentration co by solid 
lines. Respective theoretical values derived in the limit of 
high Co (see text) are shown by dashed lines. 



near the positively charged domain boundary is defined 
by the amount of carriers needed to fully compensate the 
surface bound charge and is of the order of a p /(qfCo)- 
For the concentration of Co = 1 mol% it is about 5.4 nm 
which is two orders of the magnitude smaller than the 
domain width a. This means that, at so high concen- 
trations, the electrostatic depolarization field is compen- 
sated virtually in the whole specimen volume due to 
the charge migration so that the maximum possible en- 
ergy gain and, consequently, the maximum magnitude of 
Eu, is achieved which is concentration independent. At 
smaller concentrations, especially below Co = 0.1 mol%, 
the extended space charge zone exists near the positively 
charged domain face where the depolarization field is still 
present (2lj . This makes the field suppression incom- 
plete, so that the energy gain and, consequently, En, is 



not maximum and increases with the concentration. The 
characteristic time of aging is determined, in turn, by the 
distance A the charge carriers have to cover and is ob- 
viously proportional to A/nE® ~ 1/crj. The delineated 
dependencies of both En, and r are clearly seen in Figs. 
[UandO 



IV. ANALYTICAL SOLUTION FOR MEDIUM 
AND HIGH DEFECT CONCENTRATIONS 

If the vacancy concentration is high enough in the sense 
of charge compensation mechanism discussed in the pre- 
vious section, the problem may be substantially simpli- 
fied by neglecting the diffusion contribution relevant only 
within the thin space charge zone. Then Eq. @ consid- 
ered in the drift approximation reads 

d tC = -V(/icE). (10) 

The boundary condition ((4]), where drift was out- 
weighed by diffusion, does not apply anymore. To keep 
charge balance this condition is substituted by the re- 
quirement on the surface charge density a(x,t) which 
now embraces the space charge zones and can be ob- 
tained by integration of Eq. (|10|) over an infinitesimal 
region near the boundary: 

d t a{x, t) = -q f fJ,c(x, +0, t)E z (x, +0, t). (11) 

This assumption implies that the defect concentration 
remains constant in the bulk resulting in an ansatz 

c(x, z, t) = cq + S(z)a(x,t). (12) 

Since E z (x,+Q,t) = a(x,t)/(e f +e d ), Eq-dTTJ leads to 
an equation for the surface charge density 

d t cr(x,t) = -a(x,t)/r r (13) 

with the Maxwell relaxation time r r = £n(e/ + £d)/n 
where k — q/fico is the ionic conductivity of the bulk 
material. An apparent initial condition for a(x, t) reads 
a(x,0) = <7 p sign(a/2 — x) which provides a solution 
cr(x,i) = a s {t) sign(a/2 — x) with a s (t) = a v exp (— t/r r ). 
The asymptotic condition a s (t — ► oo) = means full 
compensation of the bound charge. 

With the ansatz (fT2")l . Eq. ([3]) becomes homogeneous 
and is satisfied by the function ~E°(<r s (t)\x, z). The drift 
equation (fT0|) is then self-evident satisfied in the bulk. 
The part of the field due to charge migration equals con- 
sequently AE = -E°(ct p )[1 - exp (-t/r r )]. 

Substituting AE into Eq. ([5]) one can note that, for 
symmetry reasons, the first term in the brackets does not 
contribute to the energy . Finally, an explicit formula for 
the internal bias field follows from Eq. ©: 

t \ . 8 5 (X Pa r / i \ ^ / , \ 

E ib {t)~ -— 1 - exp (-t/7v) • (14) 

7T L EnE f 
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FIG. 4: Eit as a function of time for activation energy 
E a = 1.1 eV, acceptor concentration Co = 1 mol% and vari- 
able temperature as indicated in the plot (solid lines). En, 
averaged over the Gauss distribution of the random energies 
E a (dashed lines). 



Note that, neglecting the thickness of the space charge 
zone, the saturated (asymptotic) magnitude of En, 
achieves the concentration independent maximum value 
determined by the electrostatic energy of the stripe do- 
main structure (|7|). The aging time r is represented by 
the Maxwell relaxation time r r which is, as expected, 
proportional to 1/cq. 

The dependence (|T4|) is shown exemplary for the con- 
centration Co = 1 mol% on the Fig. [5] (dashed line). It 
describes well an increase of the corresponding computed 
curve for below the aging time but levels off at the 
magnitude approximately 15% larger than the computed 
value. We interpret this difference as a result of the grow- 
ing numerical mistake at high concentrations cq when the 
thickness of the space charge zone becomes of the order 
of (3a as discussed in Sec. IIII1 The Maxwell relaxation 
time (dashed line) describes well the aging time at all 
concentrations considered as is shown in Fig. [H 

The field En, is temperature dependent through the 
parameters e; and r r of the formula (fl"4|) . ej changes 
strongly when temperature increases towards the ferro- 
electric transition that is regarded as in Ref. Q. The 
dependence of r r is a more subtle question since, in 
the considered temperature range, the directly mea- 
sured total conductivity can be dominated by electron 
holes [13, HH Hil because of their mobility much higher 
than the ionic one. However, the density of electronic 
carriers remains, as was discussed in Sec. Ill, much 
lower than that of oxygen vacancies and cannot con- 
tribute much to the bound charge compensation. That 
is why in evaluation of r r only the ionic conductivity 
k(T) = kq exp (— E a /kT) is included with the activa- 
tion energy E a = 1.1 eV and kq — 112.8 (fi • cm) -1 for 
c = lmol% [23|. For these Co and E a , time dependences 
of En, by different temperatures are presented in Fig. 0] 



by solid lines. In fact, the activation energy for oxygen 
vacancies is not known accurately and usually estimated 
as E a = (1 ± 0.1) eV [H, Hi Assuming, as in Ref. @, 
that E a is random variable distributed with a Gaussian 
~ exp[(£; a - E) 2 /2s 2 } using E = 1 eV and s/E = 0.1, 
En, can be averaged over this distribution [H, H2] result- 
ing in a quasilogarithmic time dependencies shown by the 
dashed lines for different temperatures in Fig. 3] 

Theoretical curves in Fig. 0] represent satisfactorily ex- 
perimental time and temperature dependencies of En> 
[3, HI • Somewhat overestimated value of En, could be ad- 
justed by the factor a/ L which is, in fact, the only fitting 
parameter in this theory. In the performed calculations it 
was taken equal to 0.05 but it is, in fact, slightly depen- 
dent on the grain size and lies between 0.02 and 0.05 pa |. 
The parameter a/L may be also used to account for the 
fact that not every domain array is pinned by the local 
charges. The domain arrays which perfectly match other 
domain arrays, compensating bound charges, are stiffly 
coupled to each other. This means that for L some effec- 
tive domain length L e f f can be taken which may exceed 
L few times but is hardly larger than the grain size. 

The curves in Fig. \4\ resemble those in the theory of 
dipole reorientation [J Q though the mechanisms in- 
volved are completely different which reveals itself in dif- 
ferent dependencies on the doping concentration. Note 
that the described here space charge mechanism of do- 
main pattern fixation near domain boundaries does not 
preclude at all the dipole reorientation mechanism which 
can still be valid in the bulk of the grain. 



V. DISCUSSION AND CONCLUSIONS 

The two-dimensional model of depolarization field in- 
duced charge migration has been presented which ex- 
plains plausibly the time, temperature and doping de- 
pendencies of the internal bias field in aging ferroelecrics. 
Saturation of this field as well as descrease of the aging 
time with the increasing dopant concentration is in agree- 
ment with experimental observations @, @, • This is in 
contrast to the theory of defect dipole reorientation 0, [1] 
based on the picture of individual cage motion of oxygen 
vacancies which predicts the internal bias field propor- 
tional to and the aging time independent on the dopant 
concentration. 

A reasonable question arises: whether account of in- 
teraction between defect dipoles could modify the dipole 
rotation theory so that it explains the concentration de- 
pendence properly. This possibility should indeed be 
carefully studied. It is known, at least, that, at high va- 
cancy concentrations, substantial structural changes in 
the subsystem of vacancies may occur |44j |. This hap- 
pens, however, at about Co = 7 mol%, while the sub- 
stantial deviations from the linear dependence of E^ on 
concentration become apparent already at Co = 0.1mol% 
0, @ . Note that do pin g below 1 mol% is usually consid- 
ered as very dilute [45|]. On the other hand, one should 
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take into account possible chemical restrictions on solu- 
bility of certain acceptors in the bulk of the grains (4(| . 
In any case, a concentration of about 1 mol% is high 
enough in the sense of the presented here space charge 
mechanism of screening of polarization which allows for 
good agreement with experiment [E IS 0] • 

The proposed two-dimensional model of charge mi- 
gration can be extended to include further features and 
mechanisms which may influence aging. Note that the 
change of the 180° domain walls to the 90°-walls in the 
sketch presented in Fig. [1] does not entail strong mod- 
ification of the results and is reduced solely to the sub- 
stitution of a p = |P a | by <jp = \P s \/y/2. More impor- 
tant task is to take into account the possible change of 
the domain pattern during the charge migration. The 
metastable domain structure results from the compro- 
mise between the energy of the electrostatic depolar- 
ization field and the energy of domain walls [3J, l3t| . 
The relaxation of the electrostatic energy can trigger the 
change of the domain pattern including possible creation 
or disappearance of domain walls at the grain bound- 
aries. This process, in turn, involves the energy con- 
tribution of the mechanical stresses which has not been 
considered in this model as yet. One more idealization of 
the suggested model is an abrupt change of the polariza- 
tion at the domain boundaries. Considering space distri- 
bution of the polarization within the Ginzburg-Landau 
approach reveals that chracteristic scale at which the 
polarization gradually changes may become comparable 
with the domain width, especially when the temperature 



is not far from the temperature of the phase transition 
into the paraelectric state (47[. Including the ferroelas- 
tic interaction in the Ginzburg-Landau description may 
also substantially change the form of the domains provid- 
ing appearance of the well known needle-shaped domains 
[481 ]. These all additional features do not preclude, nev- 
ertheless, appearence of strong local depolarization fields 
which present the crucial element of the actual model of 
aging due to charged defects migration. Self-consistent 
analysis of the system evolution with many additional 
variables presents a very challenging task which may be 
addressed in the future. At the actual stage, only elec- 
trostatic arguments were observed so far which allow, 
however, comparison with the theory of defect dipole re- 
orientation where only electrostatic contributions were 
included, too. 
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